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Abstract 

Maneuvering characteristics of surface combatants in the United States Navy are often 
ignored during the design process. Key maneuvering parameters such as tactical diameter 
and turning rate are determined during sea trials after the ship enters service. In the '‘Navy 
After Next", the study of maneuvering of surface combatants will become increasingly more 
important in efforts to reduce the number of personnel required to operate the ship and 
thus reduce life cycle costs. This thesis attempts to address this issue. 

The thesis presents an Extended Kalman Filtering (EKF) algorithm to estimate the 
linear damping hydrodynamic coefficients for an Arleigh Burke Class Destroyer. Actual 
data is generated by conducting maneuvers (with a nonlinear model of the ship developed 
in a separate study) where nonlinear effects are small. The EKF then uses that data 
to estimate the hull stability coefficients (Y v . N v . Y r , and N r ) on-line in real time. The 
coefficient values determined by the EKF are then used in a simulation model and the 
results are compared to the actual trajectories. Despite the nonlinearities present in the 
actual data, the EKF provides coefficient values that reproduce trajectories with only 15% 
error. 

The linear coefficients are then used to develop simple controllers to automate maneu- 
vering for the actual ship. The parameters determined by the EKF are used to derive a 
linear time invariant (LTI) model of the ship. This LTI model then serves as the basis 
for model-based compensator designs to automatically control ship maneuvers. The first 
controller is an autopilot to regulate the ship’s heading and the second is a regulator that 
ensures the ship remains on its intended track. The performance of the compensators is 
then evaluated by simulating the performance of the LTI controllers on the nonlinear plant. 

Thesis Supervisor: Michael S. Triantafyllou 

Title: Professor of Hydrodynamics and Ocean Engineering 
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Chapter 1 



Introduction 

1.1 Background 

Iii the present naval ship design environment, reduced manning represents the cornerstone 
of any design. Maiming reductions significantly reduce the life cycle cost of the project, 
which is the main goal given today's shrinking defense budgets. Reduced manning also 
implies fewer sailors being placed in harm’s way during battle as well as enhanced quality 
of life. Reducing manning, however, affects nearly all aspects of a ship: maintenance, 
fin 'fighting and damage control, maneuvering, etc. The solution to this problem requires 
a major paradigm shift from the Navy’s current doctrine, culture, tradition, and training 
procedures. 

Naval ships must operate in a multitude of different environments and perform well in 
each. Maneuvering of the ship in each of these environments currently requires numerous 
watchstanders on the bridge and in the Combat Information Center (CIC) during normal 
operations. Of these numerous watchstanders. one must be the helmsman, the individual 
responsible for manipulating the rudder and engines to keep the ship on its desired course 
and speed. With the advances in modern control theory and computing power, the functions 
performed by many of these watchstanders may potentially be automated. This represents 
a major change in current. Navy operational doctrine, however. 

Automation of these functions requires some method of automatic control. The evolu- 
tion of control can be broken into two distinct periods. The period prior to 1957 can be 
considered the "Classical” period and the period from 1957 to the present can be considered 
tin 1 "Modern” period [33]. Classical control theory deals mainly with single-input /single- 
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output (SISO) lineal time-invariant (LTI) systems, whereas modern control theory expands 
that capability to deal with multiple-input/multiple-output (MIMO) systems. 

Automated ship maneuvering has been attempted since the invention of the gyroscope 
by E. A. Sperry in 1910. Sperry applied his gyroscope to the stabilization and steering 
of ships and later aircraft [33] using classical control methodology. Maneuvering of naval 
.surface combatants, however, is almost completely neglected in the design process. Naval 
architects are content to accept the ship's maneuvering capabilities as determined by full- 
scale sea trials after the ship has been built. The manning goal of the USN's 21 st century 
surface combatant (DD-21) is 95 personnel. Automation of maneuvering must surely be 
einploved in the design if this goal is to bo realized 

1.2 Objectives 

The main objective of this work is to apply a system identification technique known as 
the State Augmented Extended Kalman Filter (SAEKF) to identify the hydrodynamic 
coefficients of a model of the DDG-51 on-line in real time from noisy measurement data. 
The data obtained from the identification may then be used to develop simple, controllers 
that can be tested during the design process for implementation on the full-scale ship. 
It is worth noting that the simulations performed in this study may not match full-scale 
DDG-51 trial data. The model used in this study to generate the ship maneuvering data 
was developed in [43]. No attempt has been made to validate this model against full-scale 
maneuvering data. The goal is to properly identify coefficients in an assumed form of the 
maneuvering equations and attempt to reproduce the states and trajectories produced by 
the model from [43]. In terms of developing control laws, the most important factor is the 
values of the system states, rather than exact values of the hydrodynamic coefficients for a 
particular ship. In fact, linear control systems often perform quite well despite errors as large 
as 40% in the states. Therefore, if the identified coefficients produce values of the system 
states that are close to those generated by the model in [43], that form of the simulation 
equations can be assumed to be accurate and may then be used to develop control laws to 
automate the ship's maneuvering. 

In order to achieve the goal, a form of the simulation equations for ship maneuvering 
must be assumed. This form is then used to develop an SAEKF algorithm that accepts noisy 
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measurement data from ship maneuvers generated by the model in [43] and estimates the 
value of the hydrodynamic coefficients in an attempt to make the output of the simulation 
model match the actual measured data. The form of the model is based upon physical 
principles and contains the salient terms required to describe the coupled surge, sway, and 
vaw motions of the ship. Dimensionless quantities are employed throughout, the process to 
maintain generality and ensure numerical stability of the algorithm. 

The identification portion of the study limits its focus to forward motions in deep, calm 
seas with no current. The controller simulations, however, relax this condition and use 
currents and waves as disturbances to test the performance of the controller designs. 

1.3 Contributions 

This work describes and evaluates a process by which coefficients in maneuvering equa- 
tions of motion can be identified on-line in real time. From these identified parameters, 
control laws for automating the ship’s maneuvering in real time in the ship’s operational 
environment can be developed. This method can be utilized early in the design process 
on scale models to determine the effectiveness of the automatic control systems and deter- 
mine changes that need to be made. This work can be extended to apply to future naval 
combatant ship designs such as DD-21. 

This work introduces a method by which naval architects can address maneuvering 
characteristics early in the design process and develop automatic controllers to reduce the 
number of human interfaces required to maneuver the ship. The method is written en- 
tirely in the MATLAB computing environment. This work will also help to provide Ocean 
Engineering graduate students at MIT with practical experience in guidance and control 
of ocean vehicles through its potential future implementation on a scale model. Further, 
the simulations developed in the thesis may be used as the basis for future student design 
projects in ship maneuvering and control. 

1.4 Outline 

Chapter two presents a truncated derivation of the governing equations for the SAEKF 
algorithm. It then presents an application of the identification process to a very large 
crude carrier (VLCC) Esso Osaka. This is an application of the significant work conducted 
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by Professor Martin Abkowitz of MIT in the late 1970's and early 198CTs. This model 
is developed and used for verification of filter operation throughout the remainder of the 
thesis. 

Chapter three presents a brief description of the United States Navy's Arleigh Burke 
Class (DDG-51) Destroyer. It then develops the governing equations of ship motion in the 
horizontal plane. The chapter concludes by developing the assumed form of the equations 
of motion to be employed in the SAEKF. 

Chapter four develops the linear equations of ship maneuvering motion. It describes 
the results of the initial identification efforts and problems encountered. The chapter fur- 
ther describes bias and divergence in the Kalman Filter and some of the causes of these 
phenomena. The methods of addressing the problems encountered are then described. The 
chapter concludes with the results of the successful identification of the linear damping 
hydrodvnaiiiic coefficients for the DDG-51. 

Chapter five briefly describes the theory of controller design by loopshaping. Two simple 
controllers are then designed using loopshaping and the linear models produced through the 
system identification efforts. The first controller is an autopilot designed to maintain the 
ship' s actual heading about a slowly- varying reference value. Simulations of the controller 
an' performed on the nonlinear plant to determine its performance. Cross-track errors 
evidenced in the autopilot design then lead to the design of a simple track-keeping controller. 
Simulations are performed on the nonlinear plant for this design as well. 

Chapter six summarizes the work accomplished in the thesis and gives recommendations 



for future studies. 



Chapter 2 



System Identification and the State 
Augmented Extended Kalman 
Filter 



This chapter introduces the use of the Extended Kalman Filter (EKF) in system identifica- 
tion. The governing equations of the EKF appear along with a brief description of the noise 
processes. The chapter concludes with an example application of the EKF to the identifi- 
cation of unknown hydrodynamic coefficients for the very large crude carrier (VLCC) Esso 
Osaka. 



2.1 System Identification 

System Identification (SI) is the process of developing or improving a mathematical rep- 
resentation of a physical system using experimental data. This process usually assumes a 
form of a model for a physical system and adjusts the unknown parameters in that model 
to fit physical data collected from the system in question. The State Augmented Extended 
Kalman Filter (SAEKF) illustrates one technique for performing SI. The system state vec- 
tor is augmented with the system’s unknown parameters and estimated as data is collected 
from the physical system. 

Many methods exist for performing SI on mathematical models of physical systems. 
Several of these methods have been successfully applied to the ship maneuvering problem. 
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Maximum Likelihood Parameter Estimation was successfully applied in [41] to identify 
linear coefficients. Reference [13] covers several other methods such as Indirect Model 
Reference Adaptive System, Continuous Least Squares Estimation, and Recursive Least 
Squares Estimation. Abkowitz [3] successfully applied the SAEKF method to tanker ships 
in the early 1980's. The SAEKF method provides a means to update a system model on-line 
in real time. The success of the work presented in [3], coupled with the real time estimation 
capability, forms the basis for the choice of the SAEKF method in this study. 

2.2 The Extended Kalman Filter 

In the most general case, nonlinear systems and subsequent sets of discrete measurements 
of those systems can be described by a set of nonlinear, stochastic differential equations of 
the following form 1 : 



x(t) = f {x(t). t ) 4- w(t) 


(2.1) 


Ik = hk(£(<*)) + v k 


(2.2) 



where: 

t(/) represents the system state vector. 
f(£(0-f) represents the system description matrix. 

w{f) represents a zero-mean Gaussian sequence with covariance matrix Q (t), 
z k represents the measurement vector. 

hkU(/*)) represents the measurement description matrix, and 

r A . represents a zero- mean Gaussian sequence with covariance matrix R^. 

In the ship maneuvering problem, system measurements typically consist of direct mea- 
suromeiits of the state variables. Thus, the matrix h k(#(4)) reduces to a constant identity 
matrix. For example, ships often utilize a Global Positioning System (GPS) to directly 
measure the ship's position and speed. The ship also has a gyrocompass to directly mea- 
sure t lie instantaneous heading angle. The direct state measurements require no additional 
calculations to determine the value of the state. Therefore, the system measurement matrix 

l The derivations that follow appear in more detail in [19]. 



14 



reduces to a constant function equal to the identity matrix. 2 Keeping this fact in mind. 
(‘(Illation 2.2 reduces to the following: 

= h k x(t fc ) +v k (2.3) 

The equation of motion 2.1 and the measurement equation 2.3 now govern the dynamics 
of the entire system. The EKF method seeks the minimum variance estimate of x(t) as a 
function of time and the measurement data accumulated up to time t^. The sequel presents 
an abbreviated derivation of the theory behind the EKF algorithm. Reference [19] presents 
a more detailed derivation of the intermediate steps. 

2.2.1 Process and Sensor Noise 

The noise sequences w aud in equations 2.1 and 2.3 represent the process and sensor 
noise, respectively. Their properties and associated covariance matrices have yet to be 
addressed. 

In the ship maneuvering problem, the process noise represents the uncertainty in the 
assumed form of the model as well as uncertainty in predicting external disturbances. For 
example, the sea may not always exhibit a wave spectrum exactly consistent with that 
predicted by statistical data. In calm seas, however, this uncertainty is removed and only 
the uncertainty in the model form remains. Thus, the process noise covariance may be 
assumed to be constant. This is not a limitation, however, because even in a high sea state, 
the additional uncertainty due to the wave excitation will vary slowly over time and can 
be considered piecewise constant [32]. Furthermore, because the process noise represents 
uncertainty, it is assumed that the state vector and process noise are independent and 
uncorrelated random variables. 

The sensor noise represents the uncertainty in the measurements. For example, GPS 
systems available in the commercial market provide position information to within +/- 100 
feet in some instances, but not the exact location. The measurement model must address 
this uncertainty. Therefore, in the sequel, the following assumptions hold for the process 
and sensor noise: 

_> If all state variables are not measured, the matrix will consist of ones and zeros of sufficient size to 
perforin the required linear algebraic operations. 
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1. White noise processes. 3 



• I lie ship motions are slow compared to the dynamics of waves, structural vibrations, 
etc. Thus, the time constant of the process noise is much, much faster than the time 
constant of the ship motions. 

• The observation interval must be long compared to the correlation time of the sensor 
noise. Typical shipboard sensors exhibit the capability of sampling rates as high as 100 
hertz. Therefore, a sampling interval of not less than one second meets this condition 
and is more than adequate to fully capture the slow ship motion dynamics. 

2 Independent, uncorrelated, zero-mean. Gaussian random variables (denoted N(O.Q)). 



3. Covariances are constant or piecewise constant varying slowly with time. 

These assumptions are quite important in the derivations that follow. Reference [19] pro- 
vides a more detailed discussion of the noise processes. 

2.2.2 State Estimate and Error Covariance Propagation 

The EKF method seeks to minimize the error of the estimate in some statistical sense as a 
function of time. Thus, the error obviously depends on time. Define the error, x{t), and its 
associated error covariance matrix. P(t) as follows: 



where ?(/) denotes the minimum variance estimate determined by the EKF. and E[ • } 
denotes the expectation operator. From the definitions in equations 2.1, 2.4, and 2.5. it can 
be shown [19] that equations 2.6 and 2.7 govern the state estimate and its associated error 
covariance propagation. Equation 2.7 omits the time dependence of f and x for notational 

3 Xature does not exhibit white noise processes. Thus, the white noise assumption is only valid subject 
ro the conditions presented. 




= 0 for all k and t . 



x(t) = x_(t) — x(t) 

p(*j = E[x(t)m T ) 



(2.4) 



(2.5) 
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convenience. 



x(t) = f {x(t)< t) 


(2.6) 


P (t) = xp - ajf r + fx? - fx T + Q (t) 


(2.7) 



In order to solve these equations, f must be known at each instant. The value of f, however, 
depends upon the entire probability density function of x at each instant of time as shown 
bv the following relation: 

/ oc roc 

• • • / f (x. t)p(x. t) dx i . . . dx n ^ f(x, t) (2.8) 

x ./ - x 

where p(x. t) denotes the probability density function of x. Practical algorithms require 
methods of computing x and P without knowing p(x, t). One method of accomplishing this 
is by expanding f in a Taylor series about the current estimate, i. as follows (assuming the 
required derivatives exist): 

(£-$) + ■■■ (2-9) 

x=x 

Truncating equation 2.9 to include only the first order terms, taking the expected value, 
and substituting into equations 2.6 and 2.7. yields the following approximate expressions 
for the propagation of the estimated state vector and its associated error covariance: 



f(£. t) = f(x. t) + 



df_ 

dx 



x = f (ia t) 

P (t) = F(x, t)P(t) + P(t)F r (K t) + Q(*) 



( 2 . 10 ) 

( 2 . 11 ) 



where F(i\ t) is the matrix whose ij tfl 



element is defined as follows: 



, = 0MjX 

' tJ dx.j 



X=X 



( 2 . 12 ) 



2.2.3 State Estimate Update, Error Covariance Update, and the Kalman 
Gain Matrix 



The filtering equations 2.10 and 2.11 propagate the state estimate and error covariance 
over the time period between measurements < t < 1. 1 ). Therefore, to complete the 



A 
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filter and update the state estimate for the next time step, the actual state measurement 
must be taken into account. Because the EKF algorithm produces a minimum variance 
estimate of the state vector, there will likely be a difference between the propagated state 
estimate (denoted x k { — ) 4 ) and the actual measurement (denoted z k ). Thus, the updated 
state 1 estimate (denoted £*.(+)) can be computed by a linear combination of the propagated 
state estimate and the difference between the actual measurement and propagated state 
estimate'. The scaling factor in this linear combination is known as the Kalman Gain 
matrix. K^. Equation 2.13 illustrates this relation. 

■h( + ) = &(“) + K k [ 2 *. - h k (x h (-))] (2.13) 



Measurements affect the error covariance in a manner similar to the state estimate 
update. Therefore, demote the propagated error covariance as Pk( — ) and the updated 
error covariance as Pk( + ). The optimum Kalman Gain matrix, minimizes the error 
covariance update. Pk( + ). Expressing Pk( + ) as a function of of Kk yields the following 
expression to be minimized: 



Pk(+) = Pk(-) + K k £ 




h k (£*)]**(- ) r 



(2.14) 



Like the state estimate update, ^.(-b), equation 2.14 is a linear combination of the prop- 
agated error covariance and the estimation error. Minimizing equation 2.14 yields the 
following expression for the optimum Kalman Gain matrix (see [19] for details): 



Kk = -E 



£*•(-) [ia- - Ma -)] 1 



l—k - - h k (x fc )] J 



+ Rk 



-1 



(2.15) 



Note that hk depends upon the entire probability density function of x(t) similar to f in 
2.2.2.° Employing the method in 2.2.2 and expanding hk in a Taylor series about the 
current propagated state estimate. x k (— ). yields the following expression: 



Sic = hk(*jfc(-)) + H k (x k (-))(x k - x k (-)) + . . . (2.16) 



1 ( — ) and Pk( — ) represent the solutions to equations 2.10 and 2.11 on the interval th-i < t < t *.. 

’’Sot* equation 2.8. 



where: 




(2.17) 



Thus, similai to F is a matrix whose ij th element is given by the following expression: 




(2.18) 



Truncating equation 2. 1C after two terms, taking the expected value and substituting into 
equations 2.13. 2.14, and 2.15 yields the Extended Kalman Filter update equations: 



Figures 2-1 and 2-2 illustrate the process with a flowchart and timing diagram. Table 2.1 
presents a summary of the Extended Kalman Filter algorithm. Each block in the flowchart 
represents a section from Table 2.1. Figures 2-1 and 2-2 should be used in conjunction with 
Table* 2.1 to illustrate the flow of the algorithm and the governing equations for each stage 
of the process. 



M+) = -h(-) + K k \z k -h k (x k (-))] 

K k = P k (-)H k T (*,(-)) [H k (i,(-))P k (-)H k r (i,.(-)) + R k 



-1 



(2.19) 



( 2 . 20 ) 



P k (+) = [I-K k H k (f,.(-))]P k (-) 



( 2 . 21 ) 
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Figure 2-1: EKF Algorithm Flowchart 



Take a noise-corrupted 
measurement, \ 



Take a noise-corrupted 
measurement, 



- i 

x k (-) J 


{ - 

\ X k M 


Solve Plant Model ODE 
and Error Covariance ODE 


- 1 

x k + i<-> • 


i - 

i X k+] (+) 


p k<-> 


p k< + > 




p k + l« 


p k+] (+) 







l k l k+1 

Figure 2-2: EKF Timing Diagram Adapted from Figure 4.2-1 of [19] 
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Table 2.1: Summary of the Extended Kalman Filter Equations 







System Model 


x(t) — f (x(t),t) + w{t) 


Measurement Model 


Lk = h k (x(tk)) + v k 


Initial Conditions 


x(0) ~ N(xq, P o) 


Otliei Assumptions 


E — 0 V k and t 


State Estimate 
Propagation 


x = f (x.t) (Plant Simulation Model) 


Error Covariance 
Propagation 

, 


P (t) = F(x. t)P(t) + P(*)F r (i, t) + Q (t) 


State Estimate 
Update 


£*(+) = Sfc(-) + K k [z fc - h k (i fc (-))] 


Error Covariance 
Update 


Pk(+) = p-K k H k (£ fc (-))]P k (-) 


Kalman Gain 

i 

i ... . 


K k = P k (-)H k r (i,(-))S- 1 (i fc (-)) 


t 

Definitions 


F(x ( t) t) — 

a$(t) £(t)=i(t) 

H k (i*( )) = ^ . 

£=£*(-) 

S($fc(-)) = H k (x A .(-))P k (-)H k r (^.(-)) + R k 



Adapted from Table 6.1-1 of [19] 
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2.3 System Identification of the Esso Osaka 



Xo a priori information on the parameters of the DDG-51 was available in this study. Thus, 
if the identification process produced erroneous values of the coefficients, there would be 
no definitive method of isolating the problem to the model or the filter. In an effort to 
ensure proper operation of the filter prior to use in identifying the DDG-51 hydrodynamic 
coefficients, a linear model of the Esso Osaka was developed using coefficient values provided 
in [3]. 

The Esso Osaka is a 280.000 dead-weight-ton VLCC with principal characteristics listed 
m Table 2.2. Abkowitz conducted a series of experiments outlined in [3] in the late 1970\s 



Table 2.2: Esso Osaka Principle Characteristics 



LBP 


1066.3 ft 


Beam 


173.9 ft 


Draft 


92.8 ft 


A r (rudder area) 


1289.67 ft 2 


AR (rudder aspect ratio) 


1.538 


A 


314,410 Iton 


LCG (aft of midship) 


25 ft 



and early 1980's to evaluate the EKF technique as a candidate for identification of ship 
hydrodynamic derivatives. The study produced favorable results and simulations using the 
identified coefficients matched the full-scale data quite well. The successful results presented 
in [3] formed the basis for choosing the EKF method in this study. 

The mathematical calculations involved in formulating the EKF are quite involved. The 
potential for error, therefore, is rather high. Thus, in order to ensure that no mathematical 
errors had occurred, it was necessary to test the filter against a dynamical model with 
known parameters. Successful operation of the filter in identifying the known parameter 
values would indicate that the filter had been properly formulated. The values of the 
coefficients appear in Table 2.3. The negative signs have been included in the coefficient 
definitions of Y v and N v to follow the convention outlined in [49]. 
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Table 2.3: Esso Osaka Nondiiiieiisionai Hydrodynamic* Coefficients 



-y; 


0.0244 


(rn - Yr) 


0.0138 


-N v 


1.4578e-3 


(mx g - N r ) 


8.093e-4 



The results of the identification appear in Figures 2-3 and 2-4. These figures clearly 
indicate that the filter operates correctly. Notice how the parameter estimates converge to 
the* exact values* of the actual coefficients over time in Figure 2-4. This is indicative of the 
fact that the dynamical model in the filter is identical to the dynamical model producing 
the measurements. 1 lie fact that there are no modeling errors removes the requirement for 
process noise (i.e. Q = 0). Thus, the filter trusts its own state estimates and converges 
quickly to the proper values based upon the initial measurements. Figure 2-3 shows how the 
estimates of the physical states converge immediately. This is due to the exact dynamical 
model in the filter and the fact that each of the four physical states are measured. Figure 2-4 
indicates that each of the parameter estimates requires about 3 minutes to converge. This is 
because each of these states must be estimated by the filter and some dynamic information is 
required before the filter can produce good estimates. After about three minutes, the error 
between the actual and estimated state vectors has decayed to zero. Thus, the computed 
gains an* very small and the filter has effectively "learned” the dynamics of the system and 
no longer requires measurement information to produce good estimates. 
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Fi^mv 2-3: State Estimates of the Esso Osaka During a 10°/10° Zig-Zag Maneuver 




Figure 2-4: Parameter Estimates of the Esso Osaka During a 10°/10° Zig-Zag Maneim 
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Results of simulations performed using the identified coefficients appear in Figures 2- 
5 and 2-6. The identified coefficients reproduce the states and trajectory with very small 
error. This, again, can be attributed to the fact that the filter contained an exact dynamical 
model during the identification process. Note also that the trajectory shown in Figure 2-6 
is identical to the maneuver used to identify the coefficients. 
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Figure’ 2-5: Simulated States of the Esso Osaka During a 10°/10° Zig-Zag Maneuver 



Simulation of Esso OsaKa 10/10 Zig-Zag Maneuver Using Identified Coefficients 




X(yds) x10 < 

Figure 2-6: Simulated Trajectory of the Esso Osaka During a 10°/10° Zig-Zag Maneuver 
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Figure 2-7 shows the simulation results in a 10° steady turn maneuver. The simulation 
results are still quite accurate despite the differing maneuver. The error shown in the figure 
is loss than b% in the turning diameter. The error can most likely be attributed to roundoff 
errors in determining the coefficient values. This small error, however, is certainly ideal for 
control system design. 

Esso Osaka 10 Degree Steady Turn Maneuver 

5000 ; — 

4500 j- 

4000 u 

I 

3500 i- 

! 

3000 ^* 

! 

>; 2500 L 

> I 

2000 J- 

1500 L 
10001 - 
500 j- 

o 1 — 

-2000 

Figure 2-7: Simulated Trajectory of the Esso Osaka During a 10° Rudder Steady Turn 
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Chapter 3 



Nonlinear Simulation Model 



Many of the intricate details of the liullform described in the sequel have been omitted. Tins 
work uses a model developed in [43] to generate the data for the full-scale ship. Classification 
restrictions prevent verification of this model against any full-scale data for the DDG-51. 
Thus, no attempt has been made here to do so. The focus of this work lies in using the 
SAEKF process to identify the unknown parameters in the equations of motion and develop 
control systems based upon the identified parameters. The goal then becomes using the 
identified values of the parameters to reproduce ship trajectories generated by the model 
from [43]. 

The sequel develops the general equations of rigid-body motion and the assumed form of 
the simulation model to be used in the SAEKF. Note that the form of the simulation equa- 
tions developed in this study are in no way related to the form of the equations used in [43]. 
This ensures that the SAEKF does not merely estimate the parameters of a model identical 
to that which it contains. In fact, the form of the simulation equations is unimportant as 
long as the trajectory of the actual data is reproduced. Thus, the issue of coefficient can- 
cellation described in [3] becomes moot. The goal of this study is not to identify individual 
coefficient values, but to reproduce trajectories. This implies that the forces and moments 
have been adequately modeled, which is the most important factor in control system design. 

3.1 The Arleigh Burke Class Destroyer 

The United States Navy (USN) Arleigh Burke (DDG-51) Class destroyer represents the 
state of the art in operational warships today. It is a twin-screw vessel powered by four 
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LM2500 gas turbine engines coupled to controllable-reversible pitch (CRP) propellers. The 
vessel has a wider beam than conventional warships, but is still considered to have fine 
lines. Table 3.1 lists the ship's principal characteristics. The vessel uses two mechanically- 
linked rudders located downstream of the propellers to maintain course control. The ship's 
propellers maintain 100% pitch for all speeds above approximately 10 knots, and thus will 
be considered fixed for the duration of this study. The ship has port-starboard symmetry 
as do nearly all ships in existence today. 

Table 3.1: DDG-51 Principal Characteristics 



Length (LBP) 460 ft 


Beam (Bj 


59 ft 


Draft (T) 


21 ft. 


Displacement (A) 


8500 Iton 


Waterplane Area (A u ,) 


29896 ft 2 


Long Ctr Gravity (LCG) 


2.8 ft aft midship 


Block Coeff (Cfc) 


0.522 


Prismatic Coeff (C p ) 


0.615 



3.2 Governing Equations of Motion 

3.2.1 Rigid-Body Inertial Forces and Moments 

The following sections develop the equations governing the inertial forces and moments ex- 
perienced by the body in dynamic motion. Simplifications ultimately reduce the generalized 
equations to only those governing motions in the horizontal plane. 

3. 2. 1.1 Generalized Equations of Rigid-Body Motion 

The complete derivation of the generalized equations of motion for a rigid body moving on 
the surface of the earth begins by choosing an inertial coordinate system with its origin 
located at the earth's center of gravity. This is necessary to account for centrifugal forces, 
coriolis forces, etc. The body experiences these forces due to the earth’s gravitational 
field and its relative linear and angular velocity to that of the body. The effects of those 
phenomena, however, can be considered small (and thus neglected) for large objects such 
as ships. Reference [2] details this portion of the derivation, so it will not be reproduced 
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hero. 



Based upon the assumptions outlined in [2], the coordinate system is transformed into 
a body- fixed reference frame with its origin located at the midship section. Figure 3-1 
illustrates the coordinate system in a body-fixed reference frame. 

In the most general case, the 



ship may exhibit coupled mo- 
tions in all six degrees of free- 
dom shown in Figure 3-1. All 
three forces and all three mo- 
menta shown in Figure 3-1 act 
on the ship during normal op- 
erating conditions while under- 
way at sea. Equating the forces 
and moments acting on the ship 




Figure 3-1: Body-Fixed Coordinate System 



to the rates of change of linear and angular momentum experienced by the ship yields the 
generalized equations of motion shown in equations 3.1 through 3.6. 



A - m 
Y = m 
Z — rn 



u + qw - rv - x g (q 2 + r 2 ) + y g (pq - r) + z g (pr + q) 
v + rn - pw - y g (r 2 + p 2 ) + z g (qr - p) + x g (qp + r) 



■li' + pv - qv - z g (p 2 + q 2 ) + x g (rp - q) + y g (rq + p) 

I< = Ijp + (I- - I y )qr + m y g [w + pv - qu) - z g (u + ru - puo) 

M = l g q + (I x ~ Iz) r P + nl z g{u + q^ — rv) — x g {v + pv — qv) 

N — I.r + ( I y — I x )pq + m x g (v + ru — pw) — y g (w + qw — rv) 



(3.1) 

(3.2) 

(3.3) 

(3.4) 

(3.5) 

(3.6) 



where in represents the mass of the body and I x , Iy, and I z represent the moments of inertia 
about the appropriate axis. Because the longitudinal center of gravity is sufficiently close 
to the origin, and the ship is nearly symmetric, the cross-coupling inertia terms can be 
considered small and neglected. 
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3. 2. 1.2 Equations of Motion in the Horizontal Plane 



Despite the appearance of the tight coupling of all the motions in equations 3.1 through 
3.6. motions in the horizontal plane exhibit weak coupling to the out-of-plane motions [37]. 
49]. Therefore, motion in the horizontal plane can be considered separately from motions 
out-of-plane. This simplification implies the following: 

xl) = p — q — Z — K — M = 0 (3.7) 

For a ship with no list (the usual case), the transverse component of the center of gravity, 
v must lie somewhere on the x-axis. Orthogonality of the coordinate system implies the 
following: 



Vq = 0 



(3.8) 



Substituting equations 3.7 and 3.8 into equations 3.1, 3.2, and 3.6 yields the equations of 
motion in the horizontal plane: 



X = m 



u — rv — x g r 



Y — m I v 4- ru -f x g r 



N = I-r -h m 



x g (v -b ru) 



(3.9) 

(3.10) 

(3.11) 



The right-hand sides of equations 3.9 through 3.11 describe the inertial forces and mo- 
ments acting on the body during dynamic motion in the horizontal plane. The remaining 
derivations address the external forces that balance the inertial forces. 



3.2.2 Balancing Forces and Moments 

Many different forms for the equations governing the balancing forces and moments exist 
in the open literature. These models can be classified in three main categories [10]: 

1. Input-Output models. 

2. Holistic (regression) models. 

3. Modular models. 
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Input-output models describe the direct effect of varying control parameters on the maneu- 
vering response of the ship. Holistic models treat the ship as a complete entity with forces 
and moments described by a Taylor series expansion containing the pertinent kinematic and 
geometric parameters. Modular models treat all major contributing elements as separate, 
interactive modules which can be developed and tested separately. Developments over the 
last decade indicate that the flexibility of the modular approach makes it very attractive 
for future applications [10]. Thus, this type of approach has been adopted in this work. 

The following sections develop the equations governing the forces and moments 1 which 
balance the inertial forces and moments derived in 3.2.1. The total force and moment can 
hi subdivided into forces and moment* resulting from the following modular contributions: 

1. External influences (wind, waves, etc.). X ext . Y ex t . N eX f. 

2. Steady-state effects. Ao, To, Nq. 

3. Propulsive devices (propellers, thrusters, etc.), A r p, Yp, Np. 

4. Rudder forces. X p . Yp. Np. 

5. Interactions between the fluid and the hull, Xf, Yf , Nf. 

Thus, the total forces and moments may be described by the following equations: 

A = X ex t 4- AT 0 -f Ap + Xp 4- Xf 
Y = Y ext + To 4- Yp + Yp + Y f 
N — N ext ■+■ Nq 4- Np h- Np 4~ N f 

3.2.2. 1 External Forces and Moments 

This study considers maneuvering only in calm, deep water with no wind or current effects. 
This is not a likely operational environment for a full-scale ship. The goal of this study, 
however, is to reproduce maneuvering trajectories produced by the model developed in 
[43] . No environmental forces were considered in that study and therefore will not be 
considered here. Once the SAEKF method proves effective in the absence of environmental 
forces, the next logical step would be the addition of environmental forces to determine 
[ X. V. and X from equations 3.9 through 3.11. 



(3.12) 

(3.13) 

(3.14) 



31 



their effects. Therefore, X (rt . Y (T t . and N ext are taken to be zero for the duration of 
the identification portion of this work. These effects will be considered, however, during 
simulation of controller designs. This is not a limitation because a scale model may be 
tested in a pool or tank where the environment can be carefully controlled to perform the 
identification. 



3.2. 2.2 Steady-State Effects 

Steady-state effects consist offerees and moments present in steady-state motion. Thus, the 
steady-state A”o is taken to be the ship resistance at the steady-state forward speed when 
d 1 'Tier dynamic terms are zero Y and V represent steadv-state sway force and vaw 
moment, respectively. These phenomena exist primarily on single-screw ships and mani- 
fest themselves through a tendency for the ship’s stern to "walk” in a particular direction 
when the propeller thrust is small. DDG-5I has two shafts and two propellers that rotate 
in opposing directions. Thus, for this platform, the steady-state sway forces and yaw mo- 
ments cancel and can be considered zero. Furthermore, because the steady-state forces and 
moments cancel, the dynamic forces and moments will cancel as well. Thus, 

To = No = Y P = N P = 0 (3.15) 



3. 2. 2. 3 Propulsive Forces and Moments 

When a ship propels itself through water, the longitudinal force must be equal to the 
difference between the hull resistance and the propeller thrust. The open-water thrust 
provided by the propeller can be expressed as follows [49]: 



T 0 = pn 2 D 4 



K 0 + K\ J + K 2 J 2 



(3.1G) 



J (1 — w)u 
n D 



(3.17) 



where n represents the propeller speed. D represents the propeller diameter, and J represents 
the advance coefficient. The constants {Kq, K\, and K^) represent the coefficients in a 
parabolic fit of the thrust coefficient to the open- water propeller curve and (1— w) represents 
the Taylor wake fraction. Inserting equation 3.17 into equation 3.16 and carrying out the 
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requisite algebra results in the following expression: 



T (} = X P = run 2 + r) 2 nu + ry 3 n 2 



where: 



m = pD l K 0 



m = (1 - w)pD 3 K i 



V 3 = (1 - w) 2 pD 2 K 2 



(3.18) 



(3.19) 

(3.20) 

(3.21) 



Modeling of propeller thrust and torque introduces an additional state variable, n, into 
the equations of motion. Appropriate models for this variable depend heavily upon the 
type of propulsion machinery the plant contains. The DDG-51 propulsion system consists 
of four LM2500 gas turbine engines mechanically coupled to two shafts via a set of reduction 
gears. Reference [49] recommends the following form for modeling a gas turbine propulsion 
system: 



Vr^VgQE - QP 
2tvI p 



( 3 . 22 ) 



where Qe and Q p represent the engine and propeller torque, respectively, I p represents 
the polar moment of inertia for the entire engine/gear/shaft/propeller arrangement. A rep- 
resents the reduction gear ratio, and r/ r and r] g represent the propeller relative rotative 
efficiency and reduction gear transmission efficiency, respectively. This model is intended 
only to adequately model the changes in propeller speed during maneuvers. It is not in- 
tended to completely capture the performance of the engine itself. Adequate modeling of 
propeller speed changes captures the salient dynamics required for maneuvering and tra- 
jectory simulation. Therefore, this model suffices for the task at hand. Thus, the engine 
model used in this study is the same as that used in [43] 2 . 



Qe = Qi 



( afj ■ 4- b) h cf r 



n r 



Q p = prv D’ Kq(.J) = pn 2 D 5 



Q0 + Q1J + Q2J 2 



“Reference [1] provides a more detailed model to capture the performance of the engine. 



(3.23) 

(3.24) 
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whore all coefficients in equation 3.23 are taken directly from [43]. f r represents the fuel 
rate, and the constants (QorQi* and Qi) represent the coefficients in a parabolic fit, of the 
torque coefficient to the open-water propeller curve. 



3. 2. 2. 4 Rudder Forces and Moments 

Modeling the rudder forces and moments is somewhat more complicated than the derivations 
presented thus far because the angle of attack decreases due to the sway and yaw motions. 
Figure 3-2 illustrates the coordinate system and vector diagrams of the rudder forces. The 
ultimate' goal is to determine expressions for Xr and Yr. An expression for Nr then follows 

t '■< r l \ 



r 








Figure 3-2: Vector Diagram of Rudder Forces 



It is well-known that rudders are foils that act as lifting devices when positioned at some 
angle of attack. c\r . measured with respect to the oncoming flow velocity, Ur. The finite 
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I 



radius at the* leading edge of the foil requires that the lift force, L. act perpendicular to 
the oncoming flow [49]. Orthogonality of lift and drag requires that the drag force. D. act 
parallel to the flow. The lift and drag forces are defined in equations 3.25 and 3.26. 

L = -pArU%Cl{olr) (3.25) 

D = -pArUrCd{<xr) (3.26) 

When v and r equal zero, or is simply the rudder angle, 6. When the ship experiences 
sway and yaw motions, however, a/? decreases as shown in Figure 3-2 3 . The angle U 
represents the decrease in the angle of attack due to swav and vaw. Thus, the instantaneous 
lift and drag forces depend upon the instantaneous flow velocity. Ur. and instantaneous 
angle of attack, a r. The quantities shown in Figure 3-2 are determined from the relations 
shown in equations 3.27 through 3.29. 



Ur = u 2 + (v -+- x r r) 2 



(3.27) 



Or = S -f “0 



ip = arctan 



v *f x r r 



(3.28) 

(3.29) 



The lift coefficient. Cl{cvr). and drag coefficient. Cd{&r)- can be approximated with any 
appropriate model. A typical value for the drag coefficient of a foil is 0.0085 [31]. The lift 
coefficient can be modeled as a constant for small angles of attack [49]. Use of the modular 
rudder model allows the use of this approximation even in maneuvers where it appears to 
be invalid (i.e. a 20° rudder turn). This is because the actual angle of attack decreases 
rapidly with respect to the nominal as sway speed and yaw rate increase. Reference [3] used 
a 20° rudder turn to validate nonlinear coefficients identified in more violent maneuvers. 
This can be taken to mean that in that work, this type of maneuver was considered to be 
at least mildly nonlinear. The sequel will show, however, that this approximation works 
quite well for maneuvers of this type. For foils with aspect ratio ( AR ) greater than 1, [49] 

‘Figure 3-2 shows an apparent increase in the angle of attack for illustration purposes only. 0 always 
opposes to allow the forces to balance in steady state. 
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provide* the following expression for the linear lift coefficient, Cl: 



c ( 1 i 1 ■ 1 

V0.18» " r irAR ' 2n(AR) 7 

Thus, taking Cp = 0.0085 and computing Cl from equation 3.30. the lift and drag forces 
may he computed from equations 3.25 through 3.29. Once the lift and drag forces are 
known, the rudder forces are computed with equations 3.31 through 3.33 where the factor 
of two accounts for the fact that the ship has two rudders. 

Xr — 2(Ls\\\xp — Dcos'ip) 

Y r = —2 (L cos ip 4- D sin ip) 

Nr = XtYr 

3. 2. 2. 5 Hull/Fluid Interaction Forces and Moments 

The forces and moments generated by the interaction between the fluid and the hull are 
functions of several variables. Thus, Xf , Yf. and Nf may be described by the following 
equations: 



(3.31) 

(3.32) 

(3.33) 



-l 



(3.30) 



Xf = /i(u,v,r) (3.34) 

Y f = f 2 (u<v.r) (3.35) 

N f = Mu. v.r) (3.36) 

The open literature provides many expressions for f\. / 2 , and along with detailed deriva- 
tions and underlying assumptions 4 . The derivations and assumptions will not be repeated 
here with the understanding that all equations presented in the sequel conform to all as- 
sumptions described in the references from which they were taken. 

The major differences between the proposed models lie in their treatment of the non- 
linearities. Abkowitz proposed the most well-known form of the equations outlined in [37]. 
This model proposes a third-order Taylor series expansion of the fluid forces and moments 
about some steady-state forward equilibrium speed. Hwang [24] further refined this model 
based upon full-scale maneuvering trials conducted on a series of large tanker ships [3]. 
A Sov [2]. [3]. [7], [13]. [26], [52], and [54]. 
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Blanke [5 proposes a simpler form of the equations intended to capture the most important 
nonlinear it ies in terms of speed and propulsion loss. The form of the simulation equations 
chosen for this study is a hybrid of the models proposed in [5] and [24]. It should be noted 
that each of the models mentioned here are based upon the holistic method described in 
3.2.2. Because this work employs the modular method, the hydrodynamic derivatives asso- 
ciated with the rudder forces in these models 0 are captured by the rudder model described 
in 3. 2. 2. 4. Thus. fa. fa, and fa are expressed as follows 6 : 

X/ — fa ( u . v, r) = Xfrii + Xp 4- Xh (3.37) 

Y f = fa(u. v . r) = Y f ,v 4- Yf.r 4b Y h (3.38) 

Nf — fa(n,v,r) — N{,v 4- Nf-r 4- (3.39) 

Added mass and added inertia terms associated with accelerations (i.e. w, a, etc.) 
can be calculated to within sufficient accuracy for purposes of control systems through 
hydrodynamic strip theory. The terms Xf t , Y/,. and Nh in equations 3.37 through 3.39. 
however, represent some combination of unknown damping and nonlinear terms whose 
values no existing hydrodynamic theory can predict with any level of accuracy. These 
terms are currently determined through extensive model testing. Thus, the linear damping 
and pertinent nonlinear terms must be estimated to provide at least some knowledge to be 
beneficial in the design of a control system. 

3.3 Nonlinear Simulation Equations of Motion 

Completing the steps that follow results in the final form for the nonlinear simulation 
('([nations of motion. 

1. Apply the simplifications and assumptions discussed in 3.2.2. 1 through 3. 2. 2. 3. 

2. Introduce the resulting expressions for X, Y, and N to equations 3.9. 3.10. and 3.11. 

3. Solve the resulting system of equations for the state derivatives and obtain a system 
in the state-space form of equation 2.10. 

A t „>. 5 . Aa. etc. 

^Equations 3.37 through 3.39 employ standard SXAME notation for hydrodynamic derivatives. Thus, 
the coefficients in the Taylor series expansion etc. are absorbed by the hydrodynamic coefficients. 

See [37] for details. 
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Thu filial model which results from completing each of the steps listed above appears in 
Table 3.2. Coefficients of added mass and added inertia terms in Table 3.2 are determined 
from valid hydrodynamic theory and knowledge of the ship geometry. Coefficients of the 
linear damping terms and any nonlinear terms (i.e. components of A'/,, Y/ t . and 2V/, ) are to 
be treated as the unknown parameters to be estimated using the SAEKF. 



Table 3.2: Nonlinear Simulation Equations 



m — Xu 

. _ (Iz - Nr ) f 2 - (mt, - Yf)h 
0 ft 

(m- Yj,)f 3 - [rnxg - Ni)f 2 

' ~ h 

vA'IcQe - Qp 

n = 

2ixl p 

/ i = mn 2 + mnu + t] 3 u 2 + X h + X R 
h = Y h + Y R 
h = N h + Nr 

Sa — {Iz - Nf){rn - Yy) - (mx g - Yf){mx g - N ;, ) 

V\ = pD i K o 

q 2 = (1 - w)pD' i K\ 

r/3 = (1 — w) 2 pD 2 K 2 



Coefficients of inertial terms such as V r and Av to be determined from hydrodynamic theory. 
Pertinent terms in Xi,. V/, and Nh to be estimated with SAEKF. 
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Chapter 4 



DDG-51 System Identification 



Effective control system design relies upon simplified models that adequately describe the 
dominant dynamics of the plant to be controlled. These simplified models must contain the 
minimum number of parameters required to sufficiently describe the forces and moments 
acting on the plant. The equations developed in Chapter 3 are intended to describe the 
complete nonlinear plant dynamics. The dominant forces and moments are often captured 
by linear terms over a surprisingly wide range, however. Modern control system theory 
provides the capability to design controllers that are robust in the face of modeling errors. 
Thus, the first step in design is to design a controller based upon a linear model and test 
its performance on the nonlinear plant. 

The linear ship dynamics are governed by added mass and added inertia terms, as 
well as linear damping terms [37]. Hydrodynamic theory provides sufficiently accurate 
methods for calculating the former. Plant stability, however, is governed by the latter, for 
which no theory provides sufficiently accurate results. These terms must be determined 
experimentally. Several experimental methods for determining these terms are outlined in 
[37] . The drawback to these experimental methods is that they require expensive equipment 
and labor in addition to a scale model to collect the necessary data. A method of determining 
these terms using only the scale model is one of the focal points of this study. 

4.1 Linear Ship Dynamic Equations 

Ships operating at sea most often conduct maneuvers that lie within the linear regime. For 
example, ships in normal operating conditions do not usually use very large rudder angles 
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iii turns. The most common rudder angles are often between 10-15°. This type of maneuver 
lies well within the linear regime. During such maneuvers, forward speed loss and propeller 
rotational speed loss are small. They do exist, but for the purposes of control system design, 
they are sufficiently small such that they do not significantly effect the performance of the 
controller. This implies that the surge equation, u , and the propeller speed equation, n. 
in Table 3.2 are uncoupled from the sway equation, i;, and the yaw equation, f. Thus, the 
identification process reduces to the determination of and f$. 

Table 3.2 shows that and fz each contain two terms which describe the hull damping 
forces and the rudder forces. The rudder force model developed in 3. 2. 2. 4 is based solely 
upon wing theorv mid the kinematics of the problem. Thus, this force is assumed to lie 
accurately modeled and requires no identification. Therefore, the identification problem 
reduces to determining the hull damping forces, Y \ and A 7 /*. Abkowitz [37] suggests that 
m the nonlinear case, these terms are adequately modeled with a third order Taylor se- 
ries expansion. This expansion, however, contains a very large number of terms which is 
intractable from a control system perspective. For linear maneuvers, however, the expan- 
sion may be truncated after the first order terms. This truncation leads to the following 
expressions for /2 and f%\ 



ro 

It 

1 

1 

T 

1 


(4.1) 


fz = -N v v - ( mxg - N T )r 


(4.2) 



Thus, for the case of linear maneuvers, the identification problem is reduced to the deter- 
mination of the parameters —Y v , (m — y r ), —N v . and (mxg — N r ). The negative sign on Y v 
and N v was introduced to maintain the convention outlined in [49] that —Y v should have a 
large positive value. 

4.2 Initial Identification 

The 1()°/10° zig-zag maneuver proposed in [3] was used to determine the linear coefficients. 
This type of maneuver lies well within the linear regime and also provides “persistence 
of excitation" . Reference [34] states that “An open-loop experiment is informative if the 
input is persistently exciting." This effectively states that the system dynamics must be 
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continually excited in order to get any information about the parameters. A steady turn, 
for instance, would quickly reach a steady state and provide no more dynamic information 
to update the filter. Thus, persistent excitation is crucial to successful identification. 

Recall that the model generating the measurement data was taken from [43]. This model 
was adapted from a form proposed by Inoue in [26]. The model is parametrically based on 
full -scale data obtained from experiments performed over a range of operating conditions. 
The model suffers, however, during a simulated zig-zag maneuver. The implementation of 
the model in [43] suffers from a drift effect due to the rudder model. Figure 4-1 illustrates 
this effect. 
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Figure 4-1: 10"/10 w Zig-Zag Maneuver Using the Rudder Model in [43] 



Initial attempts at identifying the coefficients using this model produced biased results. 
Because no a priori information was available for the value of the unknown parameters, the 
only way to determine success or failure of identification was through simulation. Figures 
4-2 and 4-3 show the values obtained through identification and the resulting simulation. 

4.2.1 Bias and Divergence in the Extended Kalman Filter 

Divergence and bias in the Kalman Filter is well-documented in the open literature. Diver- 
gence and bias, however, do not imply filter instability. The EKF possesses some guaranteed 
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Figure 4-2: Parameter Estimates Using the Rudder Model from [43] 



Simuiatea Trajectory 




Figure 4-3: Trajectory Simulation Using the Rudder Model from [43] 
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lability properties under a few mild conditions [20]. Thus, the filter can exhibit "appar- 
ent divergence*' [12] and give incorrect (biased) parameter estimates. References [12] and 
[44] propose that bias is caused primarily by errors in the filter model. This suggestion is 
well-supported by [14] and [30]. Reference [22] provides numerical examples of the effects 
erroneous models can have on the performance of the filter. Each of these suggest that 
the model of the plant used in the filter must be sufficiently close to the actual plant to 
guarantee reliable results. 

The process noise covariance matrix, Q, used in the filter model is meant to account for 
modeling errors. Thus, the choice of Q represents a critical filter design parameter. If the 
process noise is chosen too small the filter “learns’' the wrong state values too well. This 
means that the error covariance decays too rapidly and the filter ignores any additional 
information contained in the measurements. Thus, the process noise may be said to "drive*' 
the filter and ensure that it does not place undue weighting on its own estimates. Increasing 
the process noise, however, has a strong effect on the convergence of the filter. Thus, it can 
not be chosen too high or the filter will not converge. In an effort to account for the bias, 
identifications were attempted using a wide range of process noise values ranging from 0 
and covering two orders of magnitude. Similar results were obtained in all cases. 

Much work has gone into developing an analytical method for determining the process 
noise covariance. Reference [4G] proposed a method for handling noise covariances through 
post-processing of data. References [18]. [39], and [40] all propose methods for estimation 
of the covariances in real-time. Friedland proposes another method for estimating the bias 
terms in real-time in [14]. 

4.2.2 Treatment of the Biased Estimates 

Favorable application results were demonstrated in each of the works cited above. Ap- 
plication of these techniques to the ship maneuvering problem, however, has not been 
demonstrated to the author’s knowledge. The most favorable application of the SAEKF 
technique to the ship maneuvering problem is summarized in [3], which is the culmination 
of 10 years of research conducted at the Massachusetts Institute of Technology (MIT). The 
details are outlined in the combined works of [6], [21], [24], [35], and [47]. It is evident in [24] 
that bias was a problem as well, although it is not explicitly described. In that work, the 
assumed form of the simulation model was significantly altered from the Taylor expansion 



43 



form proposed in [2] and [37] to obtain favorable results. 

This rationale has been adopted in this work as well. The drift illustrated in Figure 4-3 
most likely results from numerical problems in the model from [43] due to its parametric 
nature. The linear dynamics of ship maneuvering are well-documented and the assumed 
form shown in equations 4.1 and 4.2 are not likely to be overly erroneous. For these reasons, 
the rudder model used in [43] was replaced with the rudder model derived in 3.2. 2.4. Any 
further references to this model assume this modification has been made. 

4.3 Identification of the Linear Damping Coefficients 

4.3.1 Noise Parameters 

The magnitude of the process noise covariance matrix, Q. and the sensor noise covariance 
matrix. R. determine the weighting the filter applies to its own estimates and the measure- 
ments. respectively. During the course of the study, the trade-off between the two became 
evident. Given the sophistication of current digital modern measurement devices, the non- 
dimensional measurement noise was chosen in all cases as R = 0.01I 1 . When the process 
noise is too small, the filter “trusts" its own estimates too much and. thus, gives a very low 
weighting to the measurements (i.e. low gains). This introduces the bias discussed in 4.2.1. 
Figures 4-4 and 4-5 illustrate this phenomenon. 

The identification uses the nonlinear model from [43] to generate the data and the filter 
employs the linear model described in 4.1. Figure 4-4 shows that the states are apparently 
tracked quite closely. Note, however, that the filter never updates the surge speed, u. This 
is caused by the fact that the absence of process noise implies perfect modeling of the plant 
dynamics. Therefore, the filter has no reason to update its own propagated estimates (i.e. 
low gains). The close tracking of the sway speed, v , and the yaw rate, r, indicate that the 
nonlinear terms in these equations contribute very little during a linear maneuver. This 
indicates that the dynamics are modeled with sufficient accuracy. Thus, the filter is justified 
in trusting its own propagated estimates. 

Figure 4-5 shows the parameter estimates during the identification. Note how each of 
the parameters converges nicely to a constant value. This indicates that the error covariance 

This value is of the same order of magnitude as the non-dimensional forces, which is believed to be 
conservative. 
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Figure 4-4: State Estimates During a 10°/10° Zig-Zag Maneuver with Q=0 



Parameter Estimates with Q=0 




Figure 4-5: Parameter Estimates During a 10°/10° Zig-Zag Maneuver with Q=0 
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has decayed to zero and the filter stops updating the estimates. Recall the expression for 
the gains. K^. from Table 2.1. This equation implies the following: 

K k 0 as P k (-) -»■ 0 i k (+) = x k {~) (4.3) 



The results in Figure 4-5 are misleading, however, since they exhibit apparent divergence 
[12]. Figure 4-6 illustrates the results of the simulation after identification. The steady turn 
is chosen as the simulation maneuver to illustrate the performance of the filter estimates in 
a maneuver other than that used to identify the coefficients. Figure 4-6 shows that, while 




Figure 4-6: Simulation of a 10° Rudder Steady Turn with Q=0 



the filter does converge, the estimates produce an unstable plant. 

If perfect modeling of the system dynamics were possible, and the estimates still proved 
to be biased, the designer may introduce some “fake” noise into the system. This “fake” 
noise would continue to drive the parameters once the state estimates have converged and 
cause the filter to weight the measurements more heavily. This technique is outlined in [18] 
and suggests that a value of 0.1% of the nominal parameter value is suitable for this purpose. 
This approach fails, however, for two reasons. First, using a “nominal” parameter value 
implies some a priori knowledge of the plant. In this case, as in most cases in the physical 
world, this does not exist. Second, because the non-dimensional parameter values are small, 
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once the state estimates have converged the weighting on the measurements will be very 
small as well. Thus, the small weighting on the measurements tends to drive the filter very 
slowly. Thus, it will not converge in any reasonable period of time. Figure 4-7 illustrates 
the' effect of the “fake" noise. Note how the parameters appear to be slowly decaying with 
time. The slope is so gradual, however, that the filter will require an unacceptable amount 
of time to converge (if it does at all). 



Parameter Estimates witn Q=0 and ‘Fake" Noise 




time (sec) 



Figure 4-7: Parameter Estimates 

Given the previous considerations, the appropriate noise parameters must consist of 
some combination of process noise and “fake ' noise such that the filter continues to utilize 
the information contained in the measurements. The appropriate combination of these noise 
parameters was chosen through an iterative process of identification and simulation. Large 
process noise was required to effectively track the states. This large process noise, however, 
caused rapid convergence to the true state values. Thus, ‘Take" noise was also required to 
continue to drive the parameters to their true values. The combination of the two noise 
parameters allowed the filter to update its dynamical model, as well as obtain information 
from the measurements until a balance was achieved. The measurement noise was chosen 
as described in 4.3.1. The final values of the noise parameters are listed in Table 4.1. 
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well. The filter also does an excellent job of tracking the sway speed, v. and yaw rate. r. 
This, again, is due to the fact that sway and yaw are measured quantities and the process 
noise is high. Further, because the nonlinear effects are small in this type of maneuver, the 
dynamic model for these two terms is quite accurate. 



4.3.3 Parameter Identification 



The four additional states in the filter represent the unknown parameters to be estimated. 
Previous sections demonstrated the filter's ability to track the physical states closely, while 
giving inaccurate (biased) parameter estimates due to modeling errors and erroneous as- 
<; ri.ins ij -It# i mb.* statist i'-s. Tlu- iterative process of identification and simulation 
produced the parameter estimates illustrated in Figure 4-9. 



DDG-51 Estimated Hydrodynamic Coefficients During a 10°/10° Zig-Zag Maneuver 
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Figure 4-9: Parameter Estimates During a 10°/ 10° Zig-Zag Maneuver 



With the exception of Y v , the parameter estimates exhibit variations about some value. 
This could prove to be problematic if the goal of the identification were to establish exact 
coefficient values. The goal of control system design is to develop controllers that perform 
well with only an approximate knowledge of the coefficients. In fact, parameters that model 
the state values within 40% of the actual values often leads to acceptable control system 
performance. Fortunately, this goal is often achieved. Thus, to determine the coefficient 
value to be used in the simulation, the mean value of the last 67% of the parameter estimates 
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was takcui to be the parameter value 



4.3.4 DDG-51 Simulation After Identification 



The ultimate goal of an automatic maneuvering control system on a ship is to provide a 
reliable method of maintaining a desired trajectory. Linear control system theory has been 
proven to be quite reliable in designing control systems for plants in which nonlinearities are 
weak. Control of an inverted pendulum represents an excellent example of this assertion 2 . 
Because ship maneuvers normally lie within the linear regime, where nonlinear effects are 
small, a control system design based upon a linear ship model could very likely perform quite 
w ” m - 1-10 and 1-11 illustrate the simulation results using the identified parametei 
values. 
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Figure 4-10: Simulated States During a 10° Rudder Steady Turn 



Figure' 4-10 shows that the identified model tracks the actual nonlinear model quite 
closely in sway and yaw. It does not. however, track surge and propeller revolutions at all. 
This is because the simulation model ignores forward speed loss and propeller speed loss 
(i.e. u = n = 0). This is a good assumption because the actual speed loss in the simulation 
is about 20% and propeller speed loss is approximately 10%. These represent ideal errors 
for the application of linear theory to a nonlinear physical system. The lack of speed loss 
'See example 2.1-1 in [33]. 
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Simulated Traiectory After Identification 
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Figure 4-11: Simulated Trajectory During a 10° Rudder Steady Turn 



accounts for the higher values of sway and yaw. Higher speed in a turn implies a tighter 
turn, and thus, higher values of sway and yaw. 

Figure 4-11 shows the trajectory produced by the identified model. The error in the 
turning diameter between the identified and actual models is approximately 15%. This 
situation is ideal for the application of linear control theory to develop automatic controllers 
for the ship. Note that the identified model has a tighter turning diameter than the actual 
model. This was expected due to the higher values of sway and yaw shown in Figure 4-10. 
The identified model maintains higher speed throughout the turn. Thus, it has a smaller 
turning diameter. 
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Chapter 5 



Controller Design 



The previous chapters focused on the identification of the linear damping terms in the 
DDG-51 equations of motion. The identification process resulted in determination of the 
hydrodynamic coefficients that produce the dominant forces and moments during linear 
ship maneuvers. Thus, these terms may be employed to develop simple models which may 
form the basis for simple course- keeping controllers. 

The sequel describes the design of two controllers based upon simplified ship models. 
The first controller is a heading autopilot. The second controller regulates the lateral 
deviation (cross- track error) from the desired reference trajectory. The desired frequency 
responses are developed using a loopshaping technique. The desired control laws are then 
determined from the design loopsliape. Simulations are then performed on the nonlinear 
ship model using the linear controllers to evaluate their performance. 

5.1 Controller Design via Loopshaping 

Controller design via loopshaping is well-documented in many texts on feedback and optimal 
control system design. More detailed versions of the following presentation are described in 
[11]. [33]. and [53]. 

The overall goal of controller design via loopshaping is illustrated in Figure 5-1. The 
Loop Gain should be designed such that the closed-loop plant exhibits robust performance 
[11]. This means that the entire system (plant and controller) exhibit good performance 
and tracking at low frequencies within the system bandwidth. The system must also be 
robust to high-frequency disturbances such as sensor noise and unniodeled dynamics. The 
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Figure 5-1: The Goal of Loopshaping 



sequel describes how this process is accomplished by the control system engineer. 



5.1.1 Feedback Control 



Figure 5-2 represents a typical feedback control system with the plant disturbance input, 
(I. reflected to the plant output as in [33]. Reference [11] considers the feedback loop with 
the disturbance reflected at the plant input. The derivation of [33], however, makes the 
effect of the disturbance and desired properties of the loopshape a bit more clear. Thus, it 
is primarily that derivation that will be presented here. 
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Figure 5-2: Typical Feedback Loop 



The goal of a feedback control system is for the plant output, y, to asymptotically track 
a desired reference output, r. The plant, P, however, is acted upon by disturbances, cL 
which represent external disturbances and modeling errors. Further, the measurements of 
the output are influenced by sensor noise, n. The goal of the controller, C, is to compute 
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r!u* control input, a. that, regulates the output and maintains the desired reference value 
despite the disturbances and measurement noise. Defining the tracking error, e. as follows. 

e = r — y (5.1) 

examination of Figure 5-2 gives rise to the following expressions: 

v — r — y — n 
y = PCv — d 

0 . 1.2 System Sensitivity, Cosensitivity, and Loop Gain 

To illustrate the ideas behind loopshaping theory, one requires expressions for the closed- 
loop transfer function relations of the desired outputs, y and e. to the deterministic and 
stochastic inputs, r. d , and n respectively. Substituting equation 5.2 into equation 5.3 
results in the following expression for y 



y = PC(r — y — n) + d 


(5.4) 


(I + PC)y = PC(r - n) + d 


(5.5) 


y=(/ + PC)~' l PC{r - n) + (/ + PC)- l d 


(5.G) 


y = PC(I + PC)~' l {r - n) + (/ + PC)~ l d 


(5.7) 


where the final operation is justified because PC is square and invertible [33]. 
for the error, e. may be derived in a similar manner. 


The expression 


e = r — PCv — d. 


(5.8) 


e = r — PC{r — y — n) — d, 


(5.9) 


e = r — d — PC(e — n) 


(5.10) 


(I + PC)e = r-d + PCn 


(5.11) 


e = (7 + PC)~ l (r -d) + (I + PC)~ l PCn 


(5.12) 


e= (I + PC)~\r - d) + PC{I + PC)~'n 


(5.13) 



Equations 5.14. 5.15, and 5.16 define the system sensitivity , S', cosensitivity , T, and loop 



(5.2) 

(5.3) 
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(jam. L. respectively. 



5= (/ -rPC)~ l 


(5.14) 


T = PC{I + PC)- 1 


(5.15) 


L = PC 


(5.1G) 



Notice that 5 + T = I. Thus. S and T may not both be made arbitrarily large. Substituting 
('([nations 5.14 and 5.15 into equations 5.7 and 5.13 results in the following expressions for 
the output, y. and the tracking error, e: 

y — T(r — n) + Sd (5.17) 

e = S{r - d) 4- Tn (5.18) 

System disturbances are low-frequency phenomena such as wind and wave forces exerted 
on the ship hull. Examples of high-frequency phenomena are sensor noise and structural 
vibrational modes. The goal of the controller is to respond to low-frequency disturbances 
within the bandwidth of the system and attenuate disturbances outside the system band- 
width when* its behavior is not as well-known. For example, the controller should respond to 
wave disturbances to maintain the desired output, but should not respond to high-frequency 
structural vibrational modes whose dynamics were neglected in the system model. 

At low frequencies, equations 5.17 and 5.18 may be expressed as follows: 

y » TV 4- Sd 
e ~ S(r — d) 

At high frequencies the same equations may be represented as follows: 

y « -Tn (5.21) 

ezzTn ( 5 . 22 ) 

because the reference signal is typically a constant or at least slowly-varying. Thus, the 
system sensitivity . 5. should be small at low frequencies where the disturbances, d, are sig- 
nificant to ensure small tracking errors. Additionally, if S is small, T m 1 at low frequencies 



(5.19) 

(5.20) 
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(because 5 — 7' = 1) such that y ^ r. Conversely, the system cosensitivity. 7\ should be 
small at. high frequencies where sensor noise and other high-frequency phenomena are sig- 
nificant. If these conditions are satisfied, the controller will exhibit robust performance [ii]. 
This implies that the controller provides good tracking of reference inputs (performance) 
while simultaneously rejecting high-frequency signals where the plant model is not so well 
understood (robustness). 

5.1.3 Performance Specifications 

Frequency-dependent weighting functions are required to establish the low and high fre- 
quence performance bounds shown in Figure 5-1. These functions establish the minimum 
and maximum open-loop system gains over all frequencies. Thus, one weighting function 
should be large (i.e. large positive gains) at low frequencies to ensure good performance. 
The other weighting function should be small (i.e. negative gains) at high frequencies to 
ensure robustness. 

Determination of these weighting functions, coupled with design of the loop gain. L, 
then allows development of frequency-domain performance specifications. The idea behind 
performance specifications may be illustrated on a Nyquist plot [11] as in Figure 5-3. For 
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Figure 5-3: Nyquist Plot of Performance Specifications 



56 



good tracking. the loci of L must lie outside the disk of radius IV 1 over all frequencies. To 
ensure robustness in the face of multiplicative plant uncertainties, the disk of radius IV 2 
must not pass through the point — 1 + Oj over all frequencies. 

The goal of the controller is robust performance which implies good tracking and distur- 
bance* rejection at low frequencies and attenuation of high-frequency signals such as sensor 
noise* anel modeling errors. Recall that the system sensitivity function, 5. should be small 
at low frequencies to provide good performance and the system cosensitivitv function, T, 
shemlel be small at high frequencies to provide robustness. Considering Figure 5-3 and the 
ele'finitions e>f the quantities derived previously, it can be shown [11] that the frequency- 
domnin performance specification that ensures robust performance is given bv the following 
expression: 



\w x s\ + |iv 2 T| 



oc 



< 1 V uj 



(5.23) 



5.1.4 Design Criteria 

The length between perpendiculars (LBP) for the DDG-51 class ships is 466 feet. The main 
low-frequency phenomena in the open sea that influence the ship motions consist of waves 
and current. Currents are most often constant or at least very slowly varying. Thus, waves 
represent the phenomena that regularly contain a frequency component. 

The frequency of open-ocean gravity waves is related to their characteristic wavelength. 
The motion of the ship is also related to this wavelength. According to [49], the closed- 
loop system should respond to wavelengths of the same order as the ship LBP. Thus, the 
closed-loop system should respond to all wave disturbances of length 400 feet or greater. 
This corresponds to a characteristic circular frequency of 0.7 radians/second (Period « 9 
seconds). To allow for some flexibility in the loopshaping, the design crossover frequency, 
u;,-. should lie in the range between 0.5 - 1.5 radians/sec. 

Currents tend to "blow” the ship off course in the same way that wind blows an airplane 
off its intended course. If this disturbance is not accounted for, the ship would move far from 
its intended track. Thus, the ship must adjust its heading angle such that it accounts for the 
drift effect due to the current. Because currents are very low frequency disturbances, the 
closed-loop system must track them very accurately. Thus, the closed-loop system should 
track disturbances with an error, /?. less than 1% at frequencies below 0.01 radians/sec. 
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This corresponds to open-loop system gains of at least 40 dB below 0.01 radians/sec. Thus, 
the open-loop system gains should lie above the low-frequency performance bounds and 
below the high-frequency performance bounds illustrated in Figure 5-1 while achieving a 
reasonable closed-loop step response. The ultimate test of any controller design is how it 
operates on the actual rather than the nominal plant. The reasons for not specifying time 
domain performance criteria will be addressed in the sequel. 

5.2 Autopilot Design 

The goal of a ship autopilot is to maintain the ship's actual heading around some slowly- 
varying reference heading. The autopilot must, perform this function in the presence of 
disturbances, sensor noise, and modeling errors described in 5.1. The design speed for all 
simulations is 22 knots. This corresponds to the maximum ship speed in a plant config- 
uration with two of the four main engines on-line, one driving each shaft (split plant). 
Xoii-diniensionalization of the equations, however, allows this speed to be adjusted without 
affecting the results. The split plant configuration is required, however, to validate the 
assumption in Chapter 3 that the steady-state forces and moments due to the propeller 
cancel. 

5.2.1 Dynamic Model 

Low frequency ship motions may be described by a first-order model attributed to Nomoto 
13|. This model may be represented as follows: 

T0 + 0 = KS (5.24) 

where 0 represents the ship heading angle and represents the rudder angle. Note that 
the sway dynamics have been removed from the equation of motion. T and K represent 
the ship time constant and gain, respectively. These terms may be determined from the 
hydrodynamic coefficients identified in Chapter 4 using the procedures outlined in section 
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•‘>.3.2 of [13]. This leads to the following values foi T and K: 



T = 1.9287 (5.25) 

K = 2.04 (5.2G) 



Thus, the plant transfer function. P 0 i from input, 8, to output, ip. appears in equation 5.27 
and represents the change in heading angle to a change in the rudder angle. 




K 

s(Ts + l) 



(5.27) 



The ship's rudder contains its own dynamics. Fur instance, the rudder does not. in- 
stantaneously change its angle when commanded by the helmsman. According to [24] . the 
rudder dynamics do not play a major role in the identification process because they are 
much faster than the ship dynamics. For control purposes, however, the rudder dynam- 
ics do play a role. Classification restrictions prevent the disclosure of the actual DDG-51 
rudder swing rate. Thus, it was assumed in this work that the rudder dynamics could be 
modeled by the following transfer function. p that maps the commanded rudder angle to 
the actual rudder angle: 



p 

r sg 




(5.28) 



Tims, the new plant, P. may be represented as follows: 



p = PolP*g (5.29) 

This new plant. P. corresponds to the plant block in Figure 5-2. To guard against unrea- 
sonably large rudder angle commands during maneuvers, a rudder command angle limiter 
was added to the model as well. This limits the maximum allowable commanded rudder 
angle during maneuvers to 30°. This guards against the actual rudder angle approaching 
the physical stops on the actual ship. The goal now is to design the controller, C , using 
loopshaping to achieve robust performance. 
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5.2.2 



Loopshaping Controller Design 






111 souk' instances, it is possible to shape the loop gain, L. to obtain the desired frequency 
response. The controller may then be determined as follows: 

c = p (5.30) 

In other cases, the controller must be designed directly to obtain the desired response. The 
autopilot design followed the former procedure. The open-loop plant, P. is straight-line 
stable [37] as are many ships. Thus, the controller is not required to stabilize the plant. 

To id ve thh a physical representation, imagine the ship without any control encounters 
a current. The current will cause an initial yaw rate, but the large wetted surface area in 
the stern will develop forces and moments that counteract the initial rate of yaw. Thus, the 
yaw rate will decay back to zero after a finite period. Because the heading angle is simply 
integral of the yaw rate, or x 



'll) = r 



(5.31) 



the ship will turn to a heading that faces the current and stay on that heading when r = 0. 
Thus, the goal of the autopilot must be to ensure that the ship continues on its desired 
heading despite the current . 

Thus, the loop gain for the autopilot design may be modeled simply as a second-order 
system with an integrator whose transfer function may be represented as follows: 



— pc — UJ?i 

s{s 2 + 2(u> n + u;2) 


(5.32) 


= \/60 


(5.33) 


= 2 


(5.34) 



Figure 5-4 illustrates the results of the loopshaping design. The loopshape lias gains of 
at least 40 dB at frequencies below 0.01 radians/sec. It lies above the low-frequency per- 
formance bound indicated by W\. This indicates good tracking and disturbance rejection 
at low frequencies. The loop gains roll off at high frequencies indicating good attenua- 
tion of high-frequency signals and unmodeled fast dynamics. They further lie below the 
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Loopshape for Aufooiiot Design 




Figure 5-4: Loopshape for Autopilot Design 



high-frequency performance bound indicated by W 2 . The design also satisfies the robust 
performance requirement given by equation 5.23. 

A 11 indication of overall plant stability may be obtained from the Bode plot. Figure 5-5 
indicates a gain margin of approximately 30 dB and a phase margin of approximately 64°. 
Reference [49] indicates that a good rule of thumb for gain and phase margins are 3 dB and 
30". respectively. Thus, the autopilot design falls well within the acceptable range. 



Gain ana Phase Margins lor Autopilot Loopshape 
Gm=29 823 dB (at 7 746 rad/ sec). Pm=64 406 deg (at 0 91463 rad/sec) 




Frequency (rad/sec) 

Figure 5-5: Gain and Phase Margins for Autopilot Design 
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The frequency response of a controller is not its only measure of suitability to the task. 
The designer must also check the time domain response to ensure adequate closed-loop 
performance. Reference [36] indicates that acceptable performance for an autopilot is met 
by the following conditions: 

1. No more than 5% overshoot to a step input. 

2. 2 % settling time less than 5 ship lengths at design speed. 

Figure 5-6 indicates that the autopilot design meets these criteria. Additionally, the tracking 
error at steady-state is less than 1% as dictated by the frequency response. 



Closeo-Loop Step Response for Autopilot 

■ r, U'*i 




Figure 5-6: Closed-Loop Autopilot Step Response 



5.2.3 Closed-Loop Simulations of the Autopilot Design 

The true test of a linear controller is how well it performs in conjunction with the nonlinear 
plant from which it was derived. Incorrect assumptions on the level of uncertainty in the 
model may lead to poor performance if the unmodeled dynamics have a strong effect on the 
plant. In this case, this would arise if the nonlinearities in the equations of motion have a 
strong impact on the dynamics. These nonlinearities are believed to be weak in the ship 
maneuvering problem, however, as described in Chapter 4. This proved to be true in this 
study as well. 
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To validate the controller performance, two different simulations were developed. The 
first is a simulation of the controller's ability to reject a constant current disturbance equal 
to 5 knots. For the second simulation, a maneuver was generated by computing the required 
instantaneous heading angle required for the ship to begin on a steady course and speed 
and execute a 90° left turn and steady up on the new heading. After a period of steady 
operation, the ship then executes a 90° right turn and steadies up on its new heading. 



5.2.3. 1 Current Rejection Simulation 



Figure 5-7 illustrates the ship’s actual trajectory in relation to the dead-reckoned trajectory 
in the absence of anv disturbance. The current in this case must be simulated as an initial 
condition on the ship's sway velocity, v. due to the fact that the dynamics of position. vy. 
were neglected. The figure shows that the ship gets '‘blown 1 ' approximately 1/3 of a ship 
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Figure 5-7: Nonlinear Plant Trajectory During Current Rejection 



length off its initial track as a result of the current disturbance. It does, however, return 
to its commanded heading angle 1 , Thus, the controller takes the required action to return 
the ship to its intended heading. 

l It should be noted that the standard compass rose has been rotated 90° in each of the simulations to 
follow. Thus, a heading of 000 in the figures corresponds to a heading of 090 on a standard compass rose. 
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Figure 5-8 illustrates the change in ship heading. The ship is initially turned to the 
l ight due to the side forces and moments generated by the current interaction with the ship 
hull. The controller detects the difference between the actual and commanded headings via 
noisy measurements and generates commands to the rudder to turn the ship to the left and 
return it to its intended heading. The heading angle in the figure appears to be zero. This 
is not actually the case, however. The ship actually develops a constant sway velocity that 
results in a side-slip angle. Thus, while the figure shows the ship’s actual heading equal to 
/(to. the ship is actually turned at a small angle into the oncoming current to account for 
it. Simulations that show the actual heading angle appear in the following section where 
the dvnainics of ?/ are included in the simulation. 

Autopilot Simulation with Nonlinear Plant Model 




Figure 5-8: Nonlinear Plant Heading Change During Current Rejection 
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Figure 5-9 demonstrates the commanded and actual rudder angles during the current 
rejection simulation. The commanded rudder angles may seem excessive at first, glance. 
The figure shows, however, that during the maneuver, inclusion of the rudder dynamics in 
the equations of motion rejects these high-frequency oscillations. Thus, the actual rudder 
motion is smooth during the correction. After reaching steady-state, the rudder does exhibit 
some high-frequency oscillatory motion. The amplitude of these oscillations, however, is less 
than 0.1° over periods of several seconds. This is not deemed to be a severe condition that 
would tend to cause excessive wear on the mechanical components. 

Autopilot Simulation witn Nonunear Plant Moaei 



! Actual Ruaaer Angle 




time tsec) 

Figure 5-9: Actual and Commanded Rudder Angle During Current Rejection 

The reason for these high-frequency oscillations in the commanded rudder angle is illus- 
trated in Figure 5-10. Recall that the goal of loopshaping controller design is to ensure that 
the closed-loop combination of the plant and controller attenuate high-frequency signals. 
Because the plant has very small gains at high frequencies, the controller is free to have 
rather large gains at these frequencies. Thus, because the input to the controller is the 
noisy, high-frequency measurement of the heading angle, the controller responds in kind. 
The controller gain at the measurement frequency is approximately 8 dB. Thus, the con- 
troller actually amplifies the measurement noise. The plant gain, however, at the same 
frequency is approximately -33 dB. Thus, the effect of the noise on the entire closed-loop 
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Singula' Values of the Controller and Plant 




Figure 5-10: Plant and Controller Open-Loop Gains 



system is attenuated. 

5 . 2 . 3 . 2 Track- C hanging Maneuver 

Figure 5-11 shows the dead-reckoned and actual trajectories of the ship during a track- 
changing maneuver. The actual nonlinear plant includes the effect of forward speed loss 
due to maneuvers. Tims, because the ship loses speed during the turn and experiences 
sway motion, it is not capable of completing the tight turn based solely upon kinematic 
relations. Further, the goal of the autopilot is only to maintain the desired heading rather 
than to regulate the cross-track error. A regulator of this type is described in the following 
section. The track error is approximately three ship lengths following the initial turn. This 
is reduced, however, during the second leg of the maneuver due to the sway velocity. The 
sway velocity in the. second leg is opposite that of the initial leg. Thus, the ship actually 
moves back toward its dead-reckoned track on the second leg of the maneuver. 

Figure 5-12 shows the actual and commanded heading angle during the track-changing 
maneuver. The errors in tracking the ramp input account for the cross-track errors in the 
trajectory. The tracking error in steady-state, however, is again less than 1%. 
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Autopilot Simulation with Nonlinear Plant Model 
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Figure 5-11: Nonlinear Plant Track-Changing Maneuver with Autopilot 
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Figure 5-12: Nonlinear Plant Heading Angle during Track-Changing Maneuver with Au- 
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Figure 5-13 shows the commanded and actual rudder angles during the track-changing 
maneuver. Again, the' rudder motion is quite smooth during the major portions of the 
maneuver and is not affected by the high-frequency controller commands. Further, notice 
that the rudder sweep during the maneuver is approximately 20° over a period of about 
three minutes. Thus, the controller produces a slow, smooth rudder motion for the duration 
of the maneuver. The steady-state oscillations are again small and not deemed to cause 
excessive mechanical wear. 
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Figure 5-18: Actual and Commanded Rudder Angle During Track-Changing Maneuver 
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5.3 Cross-Track Error Controller Design 



The goal of the cross-track error controller is to regulate the ship’s lateral deviation from 
irs intended track. The previous section demonstrated that an autopilot is not sufficient for 
this purpose. The controller must meet all criteria previously stated. 

5.3.1 Dynamic Model 

To begin the design of the cross-track error controller, the coordinate system was shifted. 
In the new coordinate system, the x-coordinate represents the distance traveled along the 
reference trajectory and the y-coordinate represents the lateral deviation from the intended 
track. The heading angle now represents the deviation from the desired heading angle along 
t he intended track. The new coordinate system is illustrated in Figure 5-14. The kinematic 



relations governing ship motions in this new coordinate system is given by the following 
equations: 




X 




Figure 5-14: New Coordinate System for Track Controller 



x = u cos 'll) — v sin 0 



(5.35) 



y = u sin 0 4- v cos 0 



(5.36) 
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Linearizing and noii-diiiiensioiializing these equations leads to equation 5.37. 



V = il> (5-37) 

where the dynamics of sway have again been removed as in 5.2.1. The dynamics of x have 
also been removed since it is not a control variable. Thus, differentiating equation 5.37 
and inserting the result in equation 5.31 leads to the following expression for the control 
variable. 



V = r 



(5.38) 



Thus, using SISO loopshaping control design theory, the track controller may be designed by 
simply adding an integrator to the plant transfer function given by equation 5.29. Thus, the 
new plant transfer function that maps the control input, S c . to the lateral track deviation, 
//. becomes: 



Py = 



P 



(5.39) 



5.3.2 Loopshaping Controller Design 

Plants with double integrators can be problematic for the control systems engineer due to 
the phase shift introduced by the two poles at the origin. The two poles introduce a 180° 
phase shift in the plant at DC. Thus, because the dynamics of yaw rate, r, behave as a 
second-order system, the phase decreases as frequency increases because the plant, has slow 
polos. 

The slope of the loopshape curve as it passes through the crossover frequency (where 
tin' magnitude is one) must satisfy the condition given by equation 5.40. 

around u = u c (5.40) 

duj u r 

The phase associated with this condition is -180°. System stability requires that the phase 
he greater than this when the magnitude is one. Thus, the closed-loop system will be 
unstable if this condition is not. met. Thus, the simple loopshaping design described in 5.2.2 
may not be employed. In this instance, the controller must have some zeros to add the 
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required phase around the crossover frequency. Thus, the controller design must, proceed 
by direct design of the controller. C 2 , computation of S and T, and verification of the 
desired frequency response of the loop gain. L = P y C . Thus, the general controller transfer 
function to start the design may be defined as follows: 

c _ fc(. < j + zi)(.s + * 2 ) •••(* + *») (5 

(' s + Pl)( s + P2 ) • • • (s + pn) 

where n may be as many poles and zeros as required to achieve a desirable loopshape, 
subject to the requirement, that the transfer function remain proper (degree of numerator < 
degree of denominator) [11]. The final controller design may be described by the following 
transfer function: 



20( S + 0.5) 2 ( 5 + 6)(5 + 1000) 

(.s + ll)(.s- + 15) (.s + 30) 2 l ° j 

Figure' 5-15 illustrates the loopshape for the cross-track error controller design. The 
design proceeded iteratively by adding zeros and poles in the controller to achieve the 
desired frequency response properties. Figure 5-1G additionally shows the gain and phase 



Loopshape for TracK Controller Design 




Figure 5-15: Loopshape for Cross-Track Error Controller Design 



-’Zeros and poles may be easily added and modified using the zpk function available in MATLAB. 
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Phase (deg); Magnitude (dB) 



margins assoc iated with the design. Again, the margins fall well within the range of accepted 
practice. Notice in the phase plot of the Bode diagram how the controller adds phase to 
the system at frequencies below the crossover frequency, thus stabilizing the plant. 



Bode Diagrams 

Gm=2t.574 dB {at 7 2766 rad/sec), Pm=50 99 deg. {at 0.95135 rad/sec) 





Frequency {rad/sec; 



Finnic 5-16: Gain and Phase Margins for Cross- Track Error Controller Design 
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Figure 5-17 illustrates the system closed- loop step response. The overshoot and settling 
min' in rliis case may seem excessive when compared to conventional design practice. Coin- 
pan*. however. Figures 5-17 and 5-6. and consider equation 5.37. It would be extremely 
difficult, if not impossible, to command the dynamics of y to be faster than the dynamics 
of i' because y depends solely upon ip. Thus, its dynamics must be inherently slower. This 
physical fact is confirmed by the step response plot. Closed-loop simulation will subse- 
quently show, however, that this controller performs extremely well on the nonlinear plant 
with smooth maneuvering trajectory commands. 
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Figure 5-17: Cross- Track Error Controller Closed-Loop Step Response 
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5.3.3 Closed-Loop Simulation of the Cross-Track Error Controller 



By shifting the coordinate system to the one shown in Figure 5-14, the two simulations 
performed in 5.2.3 may be combined into one simulation. The simulation speed is again 22 
knots as in the previous simulations. The simulation steps are as follows: 

1. The ship begins on a steady course and speed with no current. 

2. After a short time, the ship experiences a current gust of 5 knots directly on its port 
beam. 

3. Shortly thereafter, the ship makes a left turn, moves two ship lengths to port and 
m readies up. 

4. After a period of steady travel, the ship executes a right turn and returns to its original 
track. 

5.3.3. 1 Maneuvering Simulation in Calm Seas with Current 

Figure 5-18 shows the desired and actual trajectories as well as the heading angle during the 
maneuver. At the beginning of the maneuver, the controller takes strong action to reject 

TraCKer Simulation with Nonimear Plant Moael 





time (sec) 

Figure 5-18: Trajectory and Heading during Nonlinear Plant Simulation 



the current disturbance. It overshoots the desired track slightly upon its return and begins 
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m compensate. The left turn then executes and the ship follows the desired track with very 
.small error until it steadies up. Notice that during the steady parts of the maneuver, the 
ship maintains a constant, heading offset to counteract the current. 

Figure 5-19 shows the commanded and actual rudder angle during the maneuver. The 
controller is again seen to respond significantly to the sensor noise for the same reason listed 
m 5.2.5. 1. The most noteworthy item in this figure, however, is the way the controller senses 

Tracker Simulation with Nonlinear Plant Model 
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Figure 5-19: Actual and Commanded Rudder Angle During Nonlinear Plant Simulation 

the current and uses it to execute the right turn to return to its original track. Note how 
the ship’s heading decreases slightly below its offset, but does not pass through zero. Thus, 
the controller uses the current to allow the ship to be “blown" back onto its desired track. 

Some may note that the ship requires a great deal of time (approximately 8 minutes) to 
complete the turning maneuvers. This may or may not be desirable. Smooth trajectories are 
required, however, to prevent the controller from over-compensating during the overshoot 
phase. This application could be extremely valuable, however, in maneuvers to avoid other 
vessels. For example, modern radar is capable of identifying contacts well before (longer 
than 8 minutes) their closest point of approach (CPA) time. Maneuvers to avoid other 
vessels are currently calculated by watchstanders using vector diagrams. Thus, a digital 
computer could calculate the CPA, compute a suitable track to avoid the contact by some 
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threshold distance, and then compute a smooth trajectory to maneuver to the new track, 
observe CPA. and maneuver back to the intended track. 



5. 3. 3. 2 Maneuver in Sea State 4 with 5-Knot Current 

The simulation would not be complete without verifying the controller’s performance in a 

sea state. The sea state is modeled as a random disturbance and filtered to preserve only the 

low-frequency components. Table G of [38] indicates that a state four sea has a significant 

wave height of 1.25-2.5 meters (4-8 feet). Additionally, Table 3 of [38] indicates that waves 

of this height represent sea conditions with the highest probability of occurrence worldwide. 

Tims, sea state four was chosen for simulation Figure 5-20 shows the wave spectrum and 

Figure 5-21 shows the time history of the wave heights. 
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Figure 5-20: Power Spectrum of Wave Model 





Figure 5-21: Time History of Wave Heights 
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The same maneuver described in 5.3.3 is performed with the current, and sea state 
introduced as disturbances. Figure 5-22 shows the actual and commanded trajectories as 
well as the heading angle during the maneuver. Notice that, similar to the maneuver in 
calm seas, the controller maintains a constant heading offset to counteract the effects of the 
current. The low-frequency wave disturbances, however, cause the ship to deviate from its 
commanded track. Thus, the ship oscillates a small amount about its commanded track. 
The amplitude of the oscillations, however, is less than 20 yards (about 15% of the ship 
lengt h). 



Tracker Simulation with Nonlinear Plant Model 
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Rudder Angle (deg) (Command (.), Actual ( -)] 



Figure 5-23 shows the commanded and actual rudder angles during the maneuver. Care- 
ful examination of Figure 5-23 indicates a low-frequency signal among the high-frequency 
measurements driving the controller. It is this low-frequency command signal that allows 
the ship to remain on track despite the disturbances. 



Tracker Simulation with Nonlinear Plant Model 




Figure 5-23: Commanded and Actual Rudder Angle during Track Change 
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Actual Rudder Angle (deg) Y (yds) 



This low-frequency component may be seen more clearly in Figure 5-24. This figure 
allows the time history during the steady portion of the maneuver between the two turns. 
Notice how the rudder corrective actions exhibit a period similar to that of the ship devia- 
tions about its intended track. In particular, the figure shows the rudder's large corrective 
action applied at the point about 1400 seconds into the simulation. This is the point where 
the ship exhibits its largest track deviation due to overshoot at the end of the turn. The 
controller senses this, however, and takes the appropriate action. 



Low-Frequency Components of Rudder Motion 
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Figure 5-24: Low-Frequency Component of Actual Rudder Angles 
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Chapter 6 



Conclusion 

6.1 Summary 

This work gives rise to a systematic: method for adding automatic maneuvering control 
system design to t lie* macroscopic ship design process. The procedure may be summarized 
as follows: 

1. Build a scale model of the ship. 

2. Outfit the model with sensors. 

.'3. Conduct a series of controlled maneuvers with the model to gather data. 

4. Design a SAEKF and use the data collected from the model to identify the ship 
hydrodynamic coefficients. 

f). Use the identified coefficients to develop controllers. 

(j. Test the controller designs on the scale model. 

This work accomplishes a piece of this process by identifying the linear damping derivatives 
and developing control laws based upon them. Specifically, this work accomplished the 
following objectives: 

• Selection of an assumed form of the equations of ship motion in the horizontal plane. 

• Design of a SAEKF to identify the linear damping hydrodynamic coefficients within these 
equations. 

• Successful identification of these coefficients which are the most critical to overall ship 
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maneuvering stability in the horizontal plane. 

• Design of two simple controllers that may be tested on scale ship models to determine 
rheii performance. 

The current method of identifying these terms involves large towing tanks and rotating 
arm basins. It also requires extensive labor to operate the equipment. These requirements 
are in addition to the need for a scale model of the ship. Thus, development of the ideas 
introduced in [3] and [43] and continued here may not only lead to an efficient, effective 
process for adding maneuvering automation to the ship design process, it may also reduce 
the infrastructure required to support it. 

6.2 Conclusions 

The following list summarizes the items deemed most noteworthy by the author over the 
course of this study. 

1. The EKF method shows definite feasibility as a method to determine linear ship 
hydrodynamic coefficients. These parameters may then form the basis for simple 
controller designs that can be tested on scale models. 

2. The* parameter estimates provided by the SAEKF need not be very accurate to pro- 
vide values suitable for controller design. The parameters determined in this study 
reproduced simulation trajectories with approximately 15% error in steady turning 
diameter. Errors as large as 40%. however, may very likely yield acceptable controller 
designs for slow maneuvers. 

3. The dynamics of sway need not be considered to develop effective linear controllers 
for the types of maneuvers considered in this study. Thus, controller designs for these 
types of maneuvers may be based on only two parameters. 

4. The EKF performs extremely well as a state estimator (due to the physical state 
measurements) despite significant inaccuracies in the filter dynamic model. Specifi- 
cally. the SAEKF model used in this study was completely linear, but still produced 
very accurate estimates of the physical states (u, u, r, and n) of the nonlinear plant 
model used to generate the measurement data. The dynamics of v and n, in fact, 
wore completely ignored in propagating the filter state estimates. 
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o. The performance of the SAEKF as a parameter estimator , however, will degrade if the 
filter model is not sufficiently accurate in describing the actual plant. The modeling 
errors will cause significant bias in the parameter estimates. 

G. The instantaneous rudder angle of attack decreases rapidly with respect to the nominal 
angle of attack as sway speed and yaw r rate increase. Thus, the linear lift coefficient 
assumption proved to be valid in the maneuvers conducted in this study despite the 
use of rudder angles traditionally believed to lie outside the linear regime. 



6.3 Recommendations for Further Study 

As in any other design study, this work leaves several areas that require further thought and 
development. Among these, the most important in the author’s mind may be summarized 
as follows: 

1. This work considers performing the identification in calm water. Environmental dis- 
turbances such as wind and waves could have a significant effect on the identification 
process. The magnitude of these effects should be investigated. 

2. Rolling motions have an effect on maneuvering in the horizontal plane. Thus, the 
inclusion of rolling motions in the nonlinear model could have an effect on the linear 
controller performance. The magnitude of these effects should be investigated. 

3. Despite the identification of four terms in this study, these four terms were reduced to 
two to develop the controllers. Identification of only these two terms using an SAEKF 
may be possible. This could make the implementation of the SAEKF much simpler. 

4. The EKF provided very accurate estimates of the surge speed and propeller rotational 
speed despite the fact that its model completely ignored their dynamics. These quan- 
tities are easily measured by the speed log and a shaft tachometer. Thus, perhaps the 
EKF could be used to provide state estimates for a MIMO controller to control the 
ship's x-coordinate as well as its y-coordinate. 

5. The proper noise covariance terms are an important design parameter in any Kalman 
filter implementation. This study determined the correct noise covariances through an 
iterative procedure of identification and simulation. Inclusion of the noise covariances 
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as filter parameters to be estimated may reduce the iterative nature of this process 
and lead to a more efficient design process. 

G. Implement the method outlined herein on the DDG-51 scale model developed at the 
MIT Towing Tank. 

7. Investigate the range over which the assumption of linearized lift on rudders is valid 
by accounting for the ship dynamics in computing the instantaneous angle of attack. 
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Appendix A 



State Augmentation of the 
Extended Kalman Filter 



The derivations presented in section 2.2 hold for estimation of the state vector describ- 
ing a nonlinear system. Augmentation of the state vector for use in estimating unknown 
parameters is the subject of the sequel. 

Recall equations 2.10 and 2.19 in 2.2.2 and 2.2.3, respectively. The variable x represents 
the estimate of the states describing the physical system. Therefore, the estimated state 
vector can be augmented with the unknown parameters in the system to allow them to be 
estimated along with the system states. In the case of DDG-51, the unkown parameters are 
constant hydrodynamic maneuvering coefficients. Maneuvering coefficients may change with 
changes in submerged geometry, speed, etc. Nondimensionalization of the equations removes 
speed considerations and submerged geometry certainly changes slowly with time. DDG- 
5 Is compensated fuel system, however, ensures that the ship maintains nearly constant 
draft over its entire range of operation. Therefore, the state vector is augmented with 
the unknown system parameters whose derivatives are zero. For example, to estimate an 
unknown parameter, h. in a system fully described by states x\ and £ 2 , the state vector and 
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Appendix A 



State Augmentation of the 
Extended Kalman Filter 



Tli<» derivations presented in section 2.2 hold for estimation of the state vector describ- 
ing a nonlinear system. Augmentation of the state vector for use in estimating unknown 
parameters is the subject of the sequel. 

Recall equations 2.10 and 2.19 in 2.2.2 and 2.2.3, respectively. The variable x represents 
the estimate of the states describing the physical system. Therefore, the estimated state 
vector can be augmented with the unknown parameters in the system to allow them to be 
estimated along with the system states. In the case of DDG-51, the unkown parameters are 
constant hydrodynamic maneuvering coefficients. Maneuvering coefficients may change with 
changes in submerged geometry, speed, etc. Nondimensionalization of the equations removes 
speed considerations and submerged geometry certainly changes slowly with time. DDG- 
51s compensated fuel system, however, ensures that the ship maintains nearly constant 
draft over its entire range of operation. Therefore, the state vector is augmented with 
the unknown system parameters whose derivatives are zero. For example, to estimate an 
unknown parameter, b. in a system fully described by states X\ and X2 , the state vector and 



84 



asM)ciate<l F matrix would be augmented as follows: 



Xi 






X 2 


— 




h 




0 



F = 



df 1 dh dh 

dx\ 3x2 83 

8 fo 8 f-2 8/2 

dx\ dx2 8$ 



000 



(A.l) 



(A. 2) 



Using this methodology, an initial estimate is made for the unknown parameters in the 
augmented stau vector and the associated error in the estimate. Tin parameters are then 
estimated in exactly the same manner as the state vector. Once the state vector has been 
augmented with the unknown parameters, and the F matrix modified appropriately, the 
EKF becomes the SAEKF. The following section presents an example application of the 
SAEKF. 
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A.l Application of the State Augmented Extended Kalman 



Filter to a Nonlinear Tracking Problem 

The Extended Kalman Filter has proven to be extremely effective for tracking applica- 
tions (such as tracking an object with a radar). The radar provides position information 
for the object, but no velocity information. Effective tracking, however, requires velocity 
information to fully describe the state of the system. 

Assume for simplicity that an object falls through space directly toward a tracking radar 
as illustrated in Figure A-l (i.e. a one-dimensional tracking problem) 1 . Assume further that 
the svstem has an unknown parameter (the ballistic coefficient, b) which affects the drag 
force on the object as it falls through space. The goal is to estimate the velocity of the 
object, v. and the unknown parameter, b , with only noisy measurements of its position, .z\ 
available as actual data. The nonlinearity arises through the dependence of drag on the 
square of the object's velocity and the ballistic coefficient ( d(b ) ex v 2 ). 



i 




Failing Body 



X 1 = X 

x2 = x = v 



X 



x3 = b 




Radar 



Figure A-l: A One-Dimensional Nonlinear Tracking Problem 



This example is taken from example G.l-2 in [19]. 
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Figure 1 A-2 illustrates the results. The upper plot shows the filter’s ability to very 
accurately estimate the object's velocity as it falls through the atmosphere. The lower plot 
illustrates the filter's ability to converge ou the true value of the ballistic coefficient using the 
measurements and computed gains. The filter does not change much initially because the 
object experiences very little drag high in the atmosphere. As the object enters the dense 
atmosphere the filter converges to the true value and remains there. Thus, one concludes 
that the filter has difficulty identifying the parameter value if it does not significantly affect 
the dynamics. 





Figure A-2: EKF Application 
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